Gyrokinetic simulations of the tearing instability 
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Linear gyrokinetic simulations covering the collisional - collisionless transitional regime of the 
tearing instability are performed. It is shown that the growth rate scaling with collisionality agrees 
well with that predicted by a two-fluid theory for a low plasma beta case in which ion kinetic dy- 
namics are negligible. Electron wave-particle interactions (Landau damping), finite Larmor radius, 
and other kinetic effects invalidate the fluid theory in the collisionless regime, in which a general 
non-poly tropic equation of state for pressure (temperature) perturbations should be considered. 
We also vary the ratio of the background ion to electron temperatures, and show that the scalings 
expected from existing calculations can be recovered, but only in the limit of very low beta. 
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I. INTRODUCTION 

The tearing instability [1] is important in magnetic fu- 
sion devices, where it drives the formation of magnetic 
islands that can significantly degrade heat and parti- 
cle confinement [2]. The related micro-tearing mode [3] 
may lead to background turbulence and is also a source 
of confinement loss. Solar flares and substorms in the 
Earth's magnetosphere are some of the many other con- 
texts where tearing plays a crucial role, inducing mag- 
netic reconnection, explosive energy release and large- 
scale reconfiguration of the magnetic field [4] . 

The evolution of the tearing instability critically de- 
pends on the relationship between the width of the cur- 
rent layer, J, where the frozen- flux condition breaks 
down and reconnection takes place, and the length scales 
characteristic of kinetic or non-magnetohydrodynamic 
(MHD) effects, such as the ion and electron skin-depths, 
di = c/ujpi and de = c/ujpe^ the ion-sound Larmor radius, 
ps = Cs/^ci^ and the ion and electron Larmor radii, pi 
and pe (see below for precise definitions of these scales; 
for Te ^ Ti^ ps ^ pi). For sufficiently large electron- ion 
collision frequency, z^g, the width of the reconnection layer 
well-exceeds all of these non-MHD scales and the mode, 
at least in the limit of strong guide-field and sufficiently 
low /3, is expected to be well described by well-known re- 
sistive MHD theory [1, 5-10]. In many plasmas of inter- 
est (e.g., the Earth's magnetosphere, modern large toka- 
maks), however, this is not the case: a decrease in the 
collisionality of the plasma leads to a decrease in the re- 
sistivity, causing the current layer width to shrink until 
it reaches or falls below the largest relevant non-MHD 
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scale; in the strong guide-field case of interest here, this 
scale is typically (for P > me/rrii) the ion-scale A^, where 
\i ^ ps ^ VPdi for /3 <C 1 and A^ ^ di for P ^ 1 and 
larger. For plasmas in which Toi/Toe = r > 1 (Toi,e are 
the ion and electron background temperatures), S < pi 
in this case, suggesting that a kinetic treatment of the 
ions dynamics is necessary. 

Given the complexity of a fully kinetic treatment, 
a variety of simplified models have been employed 
to analytically describe tearing in such cases. These 
range from cold- ion (e.^., [11, 12]) or warm- ion {e.g., 
[13, 14]) two-fluid approximations to calculations that 
include some sub-set of ion and electron kinetic effects 
(notably perpendicular ion FLR effects or electron 
Landau damping) [3, 15-18]. A further complication of 
the 6 <^ Xi regime, however, is that the tearing mode 
becomes coupled to pressure perturbations related to 
electron and/or ion diamagnetic drifts, for example, 
and most existing analytic calculations deal with this 
by invoking some type of ad-hoc closure assumption, 
e.g isothermal or adiabatic electron or ion equations of 
state. Such closures potentially play an even greater 
role as /3 approaches or exceeds unity, in which coupling 
to the sound-waves (slow waves) along the magnetic 
field is also typically important. Here we find, based 
on fully gyrokinetic simulations of the linear tearing 
mode across a range of parameters, that when two-fluid 
effects become non- negligible {S < A^), there is not, in 
general, a simple relationship between the pressure and 
density fluctuations for either the ions or electrons. The 
ratio between the two becomes a complicated function 
of position that cannot be described by simple closure 
relations. This poses a serious challenge to theoretical 
studies of collisionless or weakly collisional reconnection, 
particularly at higher /3 ^ 1 and Tgo ^ Tbi where 
S ^ pe, since a rigorous treatment would seem to require 



fully kinetic treatments of both the perpendicular and 
parallel electron and ion responses. The gyrokinetic 
ion and electron model used here allows us to explore 
numerically, over a range of plasma parameters in the 
strong guide-field limit of a simple slab geometry, the 
kinetic physics of the tearing mode and the applicability 
of some existing theories as the system transitions from 
the collisional to collisionless regimes. 



8192 (multiplied by 2/3 for dealising). Velocity space 
integrals are evaluated using Gaussian quadrature; 
the velocity grid is fixed to 20 x 16 collocation points 
in the pitch-angle and energy directions, respectively. 
Convergence tests have been performed in all runs to 
confirm the accuracy of our results. 



III. PROBLEM SETUP 



II. NUMERICAL SETUP 

We carry out simulations in doubly-periodic slab ge- 
ometry using the gyrokinetic code AstroGK [19]. The 
equilibrium magnetic field profile is 



B^oz + Bl\x)y, B,o » B^^ 



(1) 



where Bzo is the background magnetic guide field and 
By^ is the in-plane, reconnecting component, related to 
the parallel vector potential by By'^{x) = —dA^n^/dx. 
From the formal point of view, By'^ is a first-order gy- 
rokinetic perturbation. To set it up, we perturb the 
background (Maxwellian) electron distribution function 
foe{v) with a shifted Maxwellian 6fe ex v\\ {v^ v\\ are the 
velocity-space coordinates), yielding 



A^u^{x) = Amq cosh 



La:/2 



^h(x), (2) 



where A^^] = 3^3/4 (such that max|5^^(x)| = 1.0), and 
Sh{x) is a shape function to enforce periodicity [19]. The 
equilibrium scale length is denoted by a and Lx is the do- 
main length in the x-direction, set to Lx/a = 3.27r. In the 
^/-direction, we set Ly/a = 2.57r, resulting in a value of 
the tearing instability parameter A'a ~ 23.2 [10] for the 
longest wavelength mode in the system: kya = 27ra/Ly = 
0.8. Constant background temperatures (Toi^e) and den- 
sities (noi,e) are assumed for both species. We consider 
a quasi- neutral plasma, so noi = noe = ^o, and singly 
charged ions q^ = —q^ = e. 

We employ a model collision operator which satisfies 
physical requirements [20, 21] and is able to reproduce 
Spitzer resistivity [22], for which the electron- ion collision 
frequency (vq) and the resistivity {r]) are related by 



r]//j.o = 0.380uedl 



(3) 



with //o being the vacuum permeability; this formula 
holds for k±dQ < 1 {kj_^ is a characteristic scale length 
perpendicular to a mean magnetic field, kJ^^S for the 
tearing instability). Like-particle collisions are neglected. 
AstroGK employs a pseudo-spectral algorithm to 
discretize the GK equation in the spatial coordinates 
(x^y). For the linear runs reported here, it is sufficient 
to keep only one harmonic in the ^/-direction (the lowest 
harmonic is the fastest growing one). The number of 
Fourier modes in the x-direction ranges from 512 to 



We scan in collisionality and use Eq. (3) to calculate 
the plasma resistivity 7^, recast in terms of the Lundquist 
number, S = iiQaVx/r] = 2.Q'3{i'qTa)~^ {d^/ a)~'^ ^ where 
Va is the Alfven velocity corresponding to the peak value 
of By^ and ta = (i/Va is the Alfven time. Other relevant 
quantities are: 



Pi = r^/VseA/2, d, = /3-^/VseV2, 
Pe = ^'/VseA/2, 4 = /3-l/Vl/VseV2, 



(4) 



where pi^e and di^e are the ion and electron Larmor radii 
and skin-depths, respectively, a = 7X1^1 vn\^ r = Toi/Toe, 

f3e = noToe/(5^o/2Mo), PSe = CSe/^ci, ^Se = ^Tq^/tU;, 

l^ci = eBzo/^i^ Ps = PSeVTl + '^)- I^ addition to Uq, 
the adjustable parameters considered here include the 
mass ratio a, the electron beta /3e, pse/<^, and r, although 
the latter is held fixed at r = 1 except where stated 
otherwise. 

We study the collisional-collisionless transition by 
scanning in collisionality. As Ue is decreased, the different 
ion and electron kinetic scales become important. Given 
the challenge of clearly separating all the relevant spatial 
scales in a kinetic simulation, we split our study into two 
sets of runs: a smaller-pse series {pse/a = 0.02/\/2 c^ 
0.014) and a larger-pse series {pse/a = 0.2/\/2 ^ 0.14). 
Since r = 1 is held fixed except at the end of the arti- 
cle, these two sets of runs also typically correspond to 
Pi/a = 0.02 and pi/a = 0.2, respectively. 

In the former set pe, de <^ Xi ^ S <^ a; in this case 
the frozen- flux condition is broken by collisions alone, 
and since S well exceeds the collisionless electron scales 
Pe, de, such scales need not be resolved in the simula- 
tions. The ion response, on the other hand, is predom- 
inantly collisional {S > Xi) at the smallest considered 
values of 6* '^ 500 but kinetic {S < Xi) at the largest val- 
ues, S ~ 10^. Thus resistive MHD would be expected, at 
least marginally, to be valid in this case at the smaller S 
values. In the set of runs with larger- pse {pSe/d — 0.14), 
we again consider pe, de <^ Xi < a, but since pse/ct is ten 
times larger than in the previous set of runs, the ions in 
this second set are predominantly kinetic {S < A^) over 
the entire considered range of S ^ 100 — 10^. Indeed, 
at the highest values of S^ S reaches colhsionless electron 
scales {de at /3 <C 1 and pe at /3 ~ 1), and the instability 
dynamics become essentially collisionless. 

In both sets of runs, we vary S over the ranges 
mentioned above for three different sets of /3e and 



e, cr)=(0.3,0.01), (0.075,0.0025), 
These parameters are such that 
= \/T5 c^ 3.9 is held fixed and 



a = me/rrii: [(/3 
(0.01875,0.00 0625)]. 

pSe/4 = V/3e/(2a) 

thus, since ps/a is also held fixed (at either 0.014 or 
0.14), de/a is also held fixed (at either 0.0037 or 0.037, 
respectively). Given the parameter dependences of di 
and pe noted in Eq. 4, however, it is seen that the values 
of di/a and pe/a both change as Pe and a are varied 
in this manner: for pse = 0.014, di = 0.02/v^ and 
pe/ci = 0.02y^, while for pse/d = 0.14 they are ten 
times larger. 



Ca = A'a/2 (equal to 11.6 in our runs), they are 

jTAk ^ S-^^^ , 5/a - S-^^^ for /\'5/2 > 1, (5) 

-fTAk ^ 0.965;'/'cl/' , 8 /a ^ 5,"'/'ci/^ 

for A'(5/2 <C 1. (6) 

For the parameters of our simulations {Ck = 2.08 and 
Ca = 11.6) these can also be written in terms of ta and 
S as 

-fTA ^ 1.63^-^/^ , S/a - 0.78^-^/^ for A'5/2 > 1, (7) 



IV. SMALLER pse/a = 0.014 

Fig. 1 shows the tearing mode growth rate {'JTa = 
dlogA\\/dx evaluated at the X-point) and current layer 
width (full- width at half-maximum - typically somewhat 
larger than the current-profile scale-length, depending on 
the current profile) as functions of the Lundquist number 
(symbols) for T^o = T^o- Also plotted (lines) are the re- 
sults obtained from a reduced two fiuid model [12] with 
an isothermal electron equation of state. This model is 
derived under the assumption of low-/3, but exactly how 
low /3 must be for the validity of this model depends on 
how the various quantities in the model are ordered and 
is thus problem-dependent. For the ordering assumed 
in [12], it is argued that /3e <C ^/o' is required — a condi- 
tion that is marginally satisfied here only for the lowest P 
case, (/3e,cr) = (0.01875,0.000625). The two fiuid model 
is also derived under the assumption of cold ions, but 
as we show later (see fig. 4), the difference between the 
gyrokinetic results at r = and r = 1 is small. 

For the largest value of the collisionality, S < 10^, 
the gyrokinetic growth rates roll-over because the cur- 
rent layer width is too wide to satisfy the asymptotic 
scale separation, (5 <C a, assumed in the two-fiuid model 
tearing mode dispersion relation that is plotted in the 
figure. The deviation between the gyrokinetic and two- 
fiuid results at the lowest S values should therefore be 
disregarded. It is seen from the right panel that, as 
noted earlier, S > Xi for all but the largest S values 
(recall that A^, the outer- most collisionless ion-scale of 
relevance, is typically A^ ~ pse '^ pi <C (ii for /3 <C 1 and 
Ai ^ pse ^ di for Pe = 0.3). In this case, as expected, 
the two-fiuid model, at least at low-/3, recovers the well- 
known single-fiuid resistive-MHD scalings[l] and is thus 
independent of Pq. Defining the dimensionless parameter 
Ck = kyO? B' {{)) / Bq where B' {{)) is the derivative of the 
equilibrium reconnecting field at the X-point, ky is the 
linear mode- number, and Bq is the normalizing field (for 
our parameters kya = 0.8 and and aB\0)/Bo = 2.6 so 
that Ck = 2.08 for all runs), the general one-fiuid scal- 
ings are most compactly written in terms of the quantities 
TAk = ci/{VACk) = TA/Ck and Sk = CkS. Two scahngs 
are obtained depending on the product of A'^; defining 



-fTA ^ 9.145-^/^ , 6/a - 1.225-^/^ for A'6/2 <C 1. (8) 

The small and large A' expressions for 7 and 6/a both 
break-down at the point of maximum growth rate, where 
they are roughly equal: Sk ~ C^, or for our parameters, 
S ~ 600. Since S is larger or comparable to this value 
in the simulations presented here, most of our runs are 
in the either the marginally-large or small-A' regimes. 
Using for example the small-A' expressions for the value 
S = 10, 000, we obtain S/a c^ 0.03 and jta ^ 0.04, in 
rough agreement with the numerical results. The over- 
estimation of the growth rates by the two fiuid model at 
higher /3e ~ 0.3 is possibly due to either a breakdown in 
the \ow-P ordering of the fiuid model or a gradual onset 
of kinetic effects {e.g. the invalidity of a simple isothermal 
equation of state, as discussed further below). 



V. LARGER pse/a = 0.14 

We set pse/a = 0.2/>/2 = 0.14 and adjust Uq such that 
(5 < Ai, thus focusing on the regime where ion kinetic 
effects are important. The growth rate and current layer 
width versus the Lundquist number are shown in Figure 2 
(the label 6* = 00 identifies the case i^e = 0; we note that 
S may be underestimated if 6 < d^ since Eq. (3) is not 
valid in such a regime). These runs correspond to the 
same set of (/3e, a) as before; in terms of length scales, 
de/a ~ 0.037 is fixed, and di and pe change. 

As in the previous case, we observe better agreement 
between the GK and two fiuid results for lower values of 
Pe] however, as S increases and the colhsionless regime is 
approached, the agreement becomes poorer for any value 
of Pe. In this regime, electron kinetic effects (Landau 
damping and even finite electron orbits: note that for 
Pe = 0.3, S/pe ^ 2) play an important role; these are 
absent in the two fiuid model. 



VI. TEMPERATURE FLUCTUATIONS 

To better understand the discrepancies between the 
GK and two fiuid results (at high-/3e regardless of the 
collisionality and at any Pe in the collisionless regime). 
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FIG. 1: (Color online) Growth rate and current sheet width versus the Lundquist number for pse/a = 0.02/a/2. Red crosses, 
green squares, and blue circles show gyrokinetic results for (/3e, a) = (0.3, 0.01), (0.075, 0.0025), (0.01875, 0.000625), respectively. 
Red solid, green dashed and blue dot-dashed lines are the corresponding two-fluid [12] scalings. The relevant scale lengths are 
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FIG. 2: (Color online) Growth rate and current sheet width versus the Lundquist number for pse/a = 0.2/v^. Color and 
symbol schemes are the same as in Fig. 1. The scale lengths labelled on the right panel are fixed except pe (only pe for a — 0.01 
is shown). The points in the box are used to diagnose temperature fiuctuations in Fig. 3. 



we examine in Fig. 3 the validity of the isothermal clo- 
sure by diagnosing the temperature fluctuations for the 
pse/a = 0.2/a/2 and S ^7.2x 10"^ case. In Fig. 3, we plot 
the eigenfunctions of the temperature, density and elec- 
trostatic potential ((/)), and the diagonal components of 
the pressure tensor fluctuations normalized by h. These 
are measures of the polytropic indices Fg because, ac- 
cording to the polytropic law, Ps = TsToshg {s = i, e 



is a species label). We define F 



±,s 



(r. 



I yy,s)/ ^ 



and F|| s = Tzz,s^ but restrict the discussion to the par- 
allel temperature for electrons as the perpendicular com- 
ponent is at least marginally smaller (by a factor of 
0[{k±pe)^]) in Ohm's law. 

As seen in the figure, while Fy e ~ 1 outside the layer 
(thus validating the isothermal electron approximation 
in that region), it is highly peaked in the current layer 



{x ~ de) due to the Landau damping. A spatially vary- 
ing polytropic index means that the equation of state is 
not polytropic, and thus the isothermal closure is vio- 
lated. Also plotted is Fy e obtained analytically from the 
drift-kinetic electron (DKe) model [3]. We observe that 
this expression provides an excellent fit to our data for 
the case Pe = 0.01875, where electron FLR effects are 
negligible. As for the ions, Fi also varies widely over the 
ion inert ial scale in all cases, again invalidating simple 
equations of state for this species. 



VII. ION BACKGROUND TEMPERATURE 

Finally, we examine the effect of the ion background 
temperature (r), as this is the other possible source of 
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FIG. 3: (Color online) Eigenfunctions for the /3e = 0.3 case, and the polytropic indices of ions and electrons in the perpendicular 
and the parallel directions. A spatially changing polytropic index means the equation of state is not polytropic. 



discrepancy between the GK and the two fluid results. 

Plotted in Figure 4 are the growth rate and current 
layer width versus pi/pse- As before, we probe different 
values of /3e and change a such that pse/d ^ 0.14, de/a ~ 
0.037 in all cases. The Lundquist number is 5' ^ 7200. 
As in Ref. [23] , we find that the growth rate is remarkably 
insensitive to r for large P^; however, as Pe decreases we 
are able to recover the theoretically predicted scaling [16] 
of jta ~ r^/^. One possible explanation of these results 
is that, at higher /3e, coupling to sound-waves becomes 
more important and even more strongly invalidates the 
simple closure relations used in the fluid calculations. 



VIII. DISCUSSION AND CONCLUSIONS 

We have performed a set of gyrokinetic simulations of 
the linear tearing instability, varying the collisionality, 
plasma /3, mass ratio rrie/mi^ and Ti/Te. Although we 
find agreement with a two-fluid description [12] in the col- 
lisional low-/3 case, one of our main conclusions is that if 
the plasma parameters are such that two fluid effects are 
important, there is not, in general, a simple relationship 
between the pressure and density fluctuations, their ra- 
tio being a function of position. This is an intrinsically 



kinetic effect which cannot be captured by any known 
fluid closure. It seems likely that in our simulations, at 
non-small values of /3, the parallel sound wave dynamics 
become important and ion flow and temperature must be 
solved for kinetically. In addition, we find that electron 
Landau damping cannot be neglected in the collisionless 
regime. If J > pe, the effects of finite electron orbits can 
be neglected, and the Landau damping effect is analyt- 
ically tractable using the drift-kinetic model (DKe) [3]. 
For S ^ pQ the DKe model is not sufficient and a fully 
kinetic treatment is required. 

We have also shown that for /3e <C 1 (such that the 
ion sound wave and the Alfven wave are decoupled), the 
theoretically predicted dependence of the growth rate on 
r, 7 ~ r^/^ [16] is verified. However, for [3^ ^ 1 the 
growth rate in our system is remarkably insensitive to 
the background ion temperature, as noted in previous 
numerical studies [23]. 

As is widely known, linear theory breaks down when 
the magnetic island width grows beyond the width of 
the current layer, and indeed, in some physical systems, 
turbulent noise seems likely to generate seed magnetic 
islands that are nonlinear from birth. It is therefore im- 
portant to understand the importance of kinetic effects in 
the nonlinear phase — a topic not studied here. The an- 
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FIG. 4: (Color online) Ion temperature dependence of the tearing growth rate and current layer width. 



swer to this, particularly in the strong guide-field limit in 
which 3D kinetic simulations are only just starting to be 
explored, will likely depend on the aspect of reconnection 
that is of interest. As in the case of relatively small sys- 
tems without a guide field, it may be that gross features 
of the reconnection, such as the reconnection rate, can be, 
at least qualitatively, captured by fluid models. On the 
other hand, in weakly collisional plasmas, as suggested by 
the results found here, it seems likely that questions in- 
volving particle heating or energy partition, for example, 
will likely require a kinetic physics model that includes 
effects such as Landau damping and goes beyond simple 
closure schemes. The potential strength of the gyroki- 
netic model — and weakness, in some physical problems 
— is that it is designed to study the strong guide-field 
limit, in which time-resolution of the electron gyroperiod 
(necessary in particle simulations, for example, but not 
in gyrokinetics) can become a formidable challenge. Fur- 
ther work is needed to explore nonlinear reconnection in 
the strong guide- field limit, and understand the contri- 
butions that gyrokinetic simulations may offer. 
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